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Abstract 

The use of a transfer matrix method to solve the 3D Ising model is straight- 
forwardly generalized from the 2D case. We follow B.Kaufman's approach. No 
approximation is made, however the largest eigenvalue cannot be identified. This 
problem comes from the fact that we follow the choice of directions of 2-dimensional 
rotations in the direct product space of the 2D Ising model such that all eigenvalue 
equations reduce miraculously to only one equation. Other choices of directions of 
2-dimensional rotations for finding the largest eigenvalue may lose this fascinating 
feature. Comparing the series expansion of internal energy per site at the high tem- 
perature limit with the series obtained from the computer graphic method, we find 
these two series have very similar structures. A possible correct via a factor <I>(x) is 
suggested to fit the result of the graphic method. 
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1 Introduction 



Although over a half-century has passed, solving the 3D Ising model exactly is still an 
open problem. Anyone who claims to solve this model exactly should, at least, evaluate its 
partition function, internal energy per site, critical temperature and the critical exponents 
a and (3 calculated from the relevant heat capacity and magnetization per site individually. 
In addition, a crucial test, similar to one L.Onsager [ED] did in 1944, to check whether the 
results are right or wrong, is that one should compare the series expansion coefficients of 
the internal energy per site at the high temperature limit with the series obtained from 
some other methods e.g., the computer graphic method, at least up to the first three 
or four nonvanishing terms [| [0. 

Among the many various methods for deriving the partition function of the 2D Ising 
model, transfer matrix method is the oldest and original method. However, the gen- 
eralization of this method to the 3D case has had relatively little discussion. In this 
paper, we have no ambition to solve this 3D Ising model satisfying all of the requirements 
mentioned above. Instead, B.Kaufman's approach in the 3D Ising model is carried 
out step by step. Any approximation is avoided if we possibly can. In the following, it 
is shown that, when a transfer matrix formalism is set up, a spinor representation can 
work. 2-dimensional rotations in the direct product space and the feature that all of the 
eigenvalue equations reduce miraculously to only one equation also appear in the 3D Ising 
model. Even though the final high-temperature expansion series of internal energy per 
site is not exactly the same as the computer graphic method's, these two series do have 
the same structures. This discrepancy may be related to a dilemma between the choice 
of the directions of the 2-dimensional rotations in order to find the largest eigenvalue 
and losing the fascinating feature that all of the eigenvalue equations reduce to only one 
equation. Be it ever not so perfect, we hope this generalization may lay the foundations 
for further study. 

2 Transfer Matrices 

Let us consider a simple cubic lattice with I layers, each has m rows and n sites per row. 
So there are N points on the lattice, N=mnl. Periodic boundary conditions are used. 
To each lattice point, with integral coordinates r, p, (, we assign a spin variable s(r,p,£) 
which takes two values ± 1. The energy of the configuration is given by 

n m I 

E(s) = -JEE E{ s ( r > P. C)s(t+1, P , () + s(t, p, C)s(r, p+1, C)+s(r, p, C)s(r, p, C + !)}• 

T=l p=X C=l 

(1) 
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J(>0) is the coupling of a pair of neighboring spins. The partition function 

z = y.^ e{s)IT (2) 

n m I 

— E II II II e Ks ( T >p>i>)< r + 1 >p>o e Ks(T, P ,(;)s(T,p+i,o e Ks(T,p,os(T,p,(+i)^ ^ 

( s ) r=lp=lC=l 

is taken over all the 2 N possible configurations. Here K = J/T. Now we factor the 
partition function into terms each involving only two neighboring spins, giving 



Z = E ••• E <*(1, •> •) I V I 8(2, ■, ■))(s(2, , •) I V I S (3, , ■))••• (4) 

s(l, v ) s(n,v) 

(s(n,.,.)|V|s(l,.,.)> 



= TrV n , (5) 
where the matrix elements of the transfer matrix V are 

m 2 

(s(t, •, •) I V I s(r + l, •, •)) = II II e Xs(r ' p ' c)s(r+1 ' p ' c) e^ s(r ' p ''' )s(r ''' +1 ' c) e Xs(r ' p ' c)s(r ' p ' c+1) . (6) 

p=i c=i 

V can be put into a more convenient form by factoring it into the product of simpler 
matrices, 

V = V3V2V1, (7) 

m I 

(s(t, ■, •) I Vi I s(r + 1-, •)) = II II e*'(™0'( T+1 >"0 (8) 

P =i c=i 

m 2 

( S (r, , ■) I V 2 I s(r + 1-, •)) = II II S s{tMt+1iP , , (9) 

p=i c=i 

m I 

(s(r, , •) I V 3 I s(r + 1, •, •)) = II II e Ks ^Os(r, P ,c + i) 5 s{TtPt0s{T+ltPX) . (10) 

p=i c=i 

The above decomposition may be checked as follows: 

( S (t,.,.)|V 3 V 2 V 1 | S (t + 1,.,.)) (11) 

= E E E E E (s(r,;-)\V 3 \s(r+l,.,.)) 

s(r+l, v ) s'(-r, v ) s'(t+1,v) s"(r,v) s"(t+1, v ) 

(a(r + l,.,.)| S '(r,.,.)><s / (r,-,OI^|s / (r + l, V )> 
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(s'(t + 1, -,.) |/(r, •,•))(/ (r, -,.) | V 1 | s"(r + 1, •, •)) 
( s "( T+ l,.,.)| s ( r + l,.,.)) (12) 

m £ 

= IIIH E E E E E 

p=l f=l s ( T +l,p,0 s'(r,p>C) s'(r+l,p,C) s"(r,p,C) s"(r+l,p,C) 

( e ^ S(T ' P ' C)S(T ' P ' f+1) ^(r,P,C) S (r+l,p,C))(^(r+l,p,C) S '(r,p,C)) 

/ JsV(t,p,C)s'(t,p+1,C)x „ 1 \(X „ N 

l e C S '(r,p,C) S '(r+l,p,C)J l° s '(r+l, l o,C) S "(r, i0 ,C)^ 

/ JCs"(r,p,C)s"(T+l,p,C)wx U fio\ 

l e n° S "(r+l,p,C) S (r+l,p,C)^/ l 1 ^ 

m I 

— Y[ Y{ e A " s(T ' p ' c)s(T+1,p ' c) e^ s(T ' p ' c)s(T ' p+1 ' c) e i< " s(T ' p ' c)s(T ' p ' c+1) . (14) 
P =i c=i 

In the above equation, due to periodic boundary conditions, the identity 

ml ml 

n n s s(t,p,os'(t, p ,o = n n ^p+uy^p+u) 

p=i c=i p=i c=i 

is used. Furthermore, Vi, V2, V3 can be rewritten as a matrix in the direct product space. 
Observing from (8), let us define a matrix a with matrix elements 

(s(t, p, C) I a I s(r + 1, p, C)> = e ^^ +1 ^. (16) 

P K p -K 

" = < e -< ,< I (17) 

= e K (I + e- 2K a x ) (18) 



= (2 sinh2K) 1/2 e K ' a ». (19) 

7 is a 2 x 2 unit matrix. tanhK* = e~ 2K , tanhK = e~ 2K \ sintilK sinh2K* = 1. To 
simplify the matrics Vi, V2, V3, we define 

m times of I 

■ ^ — - azL 

Xij = (7 <8 7 <8 • • • <8 7) <8 • • • <8 (7 <8 ■ ■ ■ <8 o-r<8> • • • <8 7) • • • 

<g> (7 (8) 7 (8) • • • <g> 7) . (20) 

S v ' 

(=1 

Yij and Zjj are also defined similarly by replacing the Pauli matrix a x with a y and a z 
respectively. 

ml times 



V 1 = a <8 a (8 • • • <8 a (21) 
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m I 

(2sinh2K) mi/2 n n eK * x 

p=i c=i 



pX 



As for V2, we introduce another matrix b with matrix elements 

(s(r, p, I 6 I s(r + l,p,Q) = V,P,CMr + i,p,0^ s(r ' P ' CMr ' P+1 ' C) , 



(22) 



(23) 



V 2 







e 

V 

ml times 



K 






e 




A' 










\ 



-K 



6 ® 6 ' 



m i 

= n n eKZp+i xZp c 

p=i c=i 

Similarly, V3 is obtained from a matrix c, 



(24) 



(25) 
(26) 



I) (J <g> I <g) • • - J® (7*) 



(27) 



mZ times 



V 3 = 



c (X) c 1 



(28) 



= n n e ^ zp ,cZp ,c+i 

p=i c=i 

3 Spinor Representation 



(29) 



Vi, V 2 ,V 3 in 2 m '-space can be related to matrices in 2ml -spaces via Dirac V matrices. 
The process of reducing the dimensions of V had been used in the 2D Ising model. Define 
a set of matrix Y ^ satisfying anticommutation relations 



(30) 



(//,// = l,2---2m; C, C' = 1,2---Z). 



(31) 
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Every is a 2 ml x 2 ml matrix. A possible representation of is 



ri,i 


— Zi,i, 


r 2 ,i 


= Y 1A , 


r 3 ,i 


— ^1,1-^2,1 


r 4 ,i 


= -^1,1^2,1: 



(32) 
(33) 
(34) 
(35) 



r2m-l,l — -^1,1^2,1 ' ' ' X m -i t iZ m ,i, (36) 

r 2 rrt,l = -^1,1-^2,1 ••• X m _i]Y m i, (37) 

Ti,2 = ( ^1,1-^2,1 ••• X mj i, ) Zi t 2 = U±Zi t 2, (38) 

r 2 , 2 = (X ltl X 27l ---X m7l )Y lt2 = U 1 Y 1 , 2 , (39) 

r liC = u 1 u 2 ---u c - 1 z ltC , (40) 

r 2 , f = U 1 U 2 ---U C - 1 Y 1X , (41) 

r 2m _i,/ = U 1 U 2 ---U l - 1 Z miU (42) 

r 2m ,i = UJJi- ■ -Ui-iY mh (43) 

where 

U c = (I <g) • • ■ <S> /) <S> • • • <8> ( <J X g cr x ® ■ ■ ■ Q ® (/ <8> ■ • • <8> I) • • ■ <8> (/ <8> • • • <S> I). (44) 

( th block 

The total number of the T matrix is 2ml. A special 2 ml x 2 ml matrix U is defined as 

U = & x <g> o x <S> • • ■ <8> o x (45) 

S v ' 

ml times 
I 2m 



= i ml 



n(n^c)- (46) 

C=l ix=l 



U and Uc have the following relations: 



U 2 = I, U 2 C = I, (47) 

U(I + U) = I + U, U C {I + Uc) = I + U c , (48) 

U(I -U) = U -I, U C (I -U c ) = U C - 1, (49) 

{U,T^} = 0, {U o T^} = 0, (50) 
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By definition of T, we notice that 



then 



Simlarly, 



then we have 



r 2 p,cr 2 p-i, c - y pX z pX - ix P £ , 



Vi = (2sinh2K) ml/2 [] IT e iK ' r *»-™ r ™ 
p=i c=i 



i m— 1 

y 2 — Y]_{e KZm ' (Zl '' : Y\ e KZp ' <:Zp+1 ' c } 
C=i p=i 

£ m— 1 

= fI{ e iia/ < ri -< r2m -< JJ e ^r 2p , c r 2p+1 , c ^ 
C=i p=i 



With the identity, 



(51) 

(52) 
(53) 
(54) 

(55) 
(56) 

(57) 



V2 can be rewritten as 



m— 1 



V 2 = n{[-( J + f/c) e ^ ri,cr2m ' c + -( J - U c )e~ iKr ^ r ^] [] e^ 2 "^^ 1 ^}. (58) 
C=i 2 2 p=i 



Since commutes with I^p^I^p+i^, the projection operators \ (I ± £7^) project V 2 into 
2* pieces. For 

m l—l 

V 3 = ft (e KZ »^ J] e^^ 1 ) , (59) 
the situation seems more complicated. 



Z p,C Z p,(+i 



( Z p,c Z p+i,d( Z p+i,<: Z P+u) ■ ■ ■ { z m,c, z u+i){ z u+i z 2,c+i) 

■•■( z p -i,t+i z p,c+i), (6°) 

( - i^2p+l,C^2p,c)( - ^r2p+3,C r 2p+2,c) • • • 

( - i^i,c+i^2m,c) • • • ( - ir 2 p-i^ + ir 2 p-2,c+i), (6i) 

^ m r2p^r2p+i^r2 P+ 2,c • • • ^2m,c^i,c+^ ■ ■ ■ r2p-i,c+i, (62) 

i m ^2p,c{^2 P +i^ ■ ■ ■ r 2 p-2,c+i)r 2 p-i,c+i, (63) 

+iT 2 p,(W2p+i^2p-i,(,+i-, (64) 

(65) 
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W 2 p+i£ is defined as 

WVh,c = i m ~ 1 T2p+i,(F2p + 2,<; ■ ■ •r 2 p-3,c+ir 2 p-2,c+i (66) 

(p+W (p-i,C+i) 
= J® •••/(g) (g---<g / o^ N <g/---<g/. (67) 

V v ' 

m— 1 times of cr x 

^2p+i,c has the property that it anticommutes with inside the region that the integral 
coordinates (/i, a) from (2p + 1,C) to ( 2p — 2,£ + 1 ), whereas it commutes with r MjQ , 
outside of that region. 

z p,i z p,i = (^^+i,«) ••• ( Z m,i z i,i)( z i,i Z 2,i) ■■ -(Zp-x^Zp^) (68) 

= — i m UT 2 p t iT 2p+ i : i ■ ■ ■ r 2mi ir ljl r 2) i • • ■ r 2p _ ljl (69) 

= +iUW2p+i : iT2p-i t iT 2 p : i- (70) 

W2 P +i,i = i m 1 r 2 p + i^r 2 p +2i i • • • T2m,i^i,i • • • r 2p _ 2j i (71) 

= a x <g • • • <g (7s <gf (g • • • <g / (g <g • • • (g a x . (72) 
(1,1) (p-1,1) (p+i,0 (m,o 

Then we have 

m, i— 1 

y 3 = JJ| e +^^2 P+ i,;r 2p -i,ir 2p ,, -q e «w H1 , c r 2p , { r 2 ,_ u+1 ^ ^g-j 
p=i C=i 

A remarkable observation of Kaufman is that decomposing V of the 2D Ising model 
into the product of factors like e^ rr , which is interpreted as a two-dimensional rotation 
with rotation angle 6 in the direct product space. We follow this spirit and decompose 
the factors into several 2D rotations, 

e +iKUW 2p+1;l r 2p _ 1A r 2p;l 

— I(J _|_ U^ e +i-KW 2 p + i,iT2p-i,iT 2 p,i _|_ l^j _ jj-j e -iKW 2p+ i : iT 2p -i,iT 2pil 

= \{I + U){ i(J + W 2p+ i,i)e +lKr *»-^ r ^ + W 2 p + i, l )e- lKr ^ r ^ } 

~ U) + W 2p+1 , l )e- iKr ^ r ^ + 1(7 - ^p+i^e^ 2 "- 1 ^ 2 ^ } (75) 



e 



*-^W / 2p+l,C r 2p,C r 2p-l,C+l 

= i(J + ^ 2 p + i,c)e lXr2 "^ r2 "-^ +1 + i(J - ^ 2 p+i, c )e- iXr2 -c r2 -^+ 1 . (76) 

In the 2D Ising model the V are decomposed into 2 pieces. This is not so simple in the 
3D Ising model. Due to the projection operators \ (I ± U), \ (I ± U^), |(I ± W), V has 
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2 l x (2 x 2 ml ) pieces. Only one piece will produce the largest eigenvalue, which dominates 
the value of partition function. 

Since W^p+i^ ma y commute or anticommute with T^ a , we have to check the com- 
mutation relations between the projection operators and the product of the factors like 
e^ err , since the dogma states that two matrices are simultaneously diagonized if and only 
if these two matrices commute. V comprises the product of projection operators and a 
lot of e5 err . So if the commutation relations are not valid, the whole scheme may break 
down. Commuting |(7 ± U) with the product of all factors like e^ err does not give any 
trouble. Any single projection operator |(7±C/^) or |(J± W) may do not commute with 
some e^ 9rr . However, fortunately, the product of all the projections in any one piece of 
V, e.g. 

i l i mii 

p = 9 a + u)( n + u,)) (nn^+ M- ^ 

z C'=i p=ic=i z 

do commute with any e5 9rr , because e5 err passes through P will change the sign of ^OTT 
even number times and it does not change sign eventually. 

How to choose the proper piece and obtain the largest eigenvalue? We have no an- 
swer. A possible rule may be followed that all eigenvalue equations should reduce to one 
equation. We will show this point in the next section. This is a beautiful feature in the 
2D Ising model. Let us consider one possible piece of V, 

V = (2sinh2K) ml/2 V' , (78) 

m I— I 

p=i c=i 

l m—l 
{ J| e iKr i,< ,T 2m,c' 7J e ^ r 2p,C r 2p+l,C} 

C=i P '=i 

l m 

{ n II e tfc *V-i,c" r 2p",c"}. (79) 
C"=ip"=i 

In essence, V includes the repetition / times of the same rotations as in the 2D Ising 
model, appearing in the second and third brackets of (79), and the new rotations in the 
first bracket of (79), relating to the third dimensional coupling beyond the 2D Ising model. 

4 Eigenvalue Equations 

The rotation operator in the spinor representation, 

Sx,(0) = e^ r * (A ^ a), (80) 
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has a one-to-one correspondence to the 2D rotational matrix oj(\<j \ 9) of the V matrix. 

s^(9)r a s Xa (9) = 5>(a<t 1 0)™^, (si) 

K 

S^(9)T x S Xa (9) = r A cos^ + r a sin^, (A ^ a), (82) 

S^(9)r a S Xa (9) = -T x sm9 + r a cos9, (a ^ A), (83) 

S^(9)r a S Xa (9) = F a , ( a ^\a^a). (84) 

S Xcr (9) is a 2D rotations in the direct product space. The rotations, e ± ^ rr , have different 
rotational angles, +9 and —9. We mean they have different directions of rotations. The 
eigenvalues of the rotational matrix u are 1 with (2ml — 2)-fold degeneracies and e ±ld two 
nondegenerate eigenvalues, whereas the eigenvalues of S Xa in the spinor representation 
are e ±l ^ each with 2 mi_1 -fold degeneracies. 
The correspondence with V is 



(2sinh2K) 



ml/2 



(85) 



u>' = {Y[u(2p-l,l;2p,l\+2iK)l[u(2p,{;2p-l,{ + l\2iK)} 
p=i c=i 

I m— 1 

{ n ^(i, C'; 2m, c' I 2^) J] w(2p, C; 2p + 1, C I 2iK)} 

C'=i p'=i 

I m 

iU II ^(2p"-l,C";2p",C"|2^*)} (86) 

a" 1 " 1 
C =1/9 =1 



(87) 



( n 



V 

fit 




fi 

fit 



••• 

-fi 





fi 




fi 



fit A 
■■■ fit 



n 



J 



2mlx2ml 



-fit\ 





fi 



(88) 



(89) 



2rax2ra 
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A = 



fit 



CC IS c 

-is*c c*c 



B 



i 

2 

1+c* 



«s 



_ 1 

(4^) 



(90) 



(91) 



s = sinh 2X, c = cosh 2K, s* = sinh 2K\ c* = sinh 2K* . (92) 

The matrix Q is just the same matrix considered in the 2D Ising model. ujIuj^ujI is in 
symmetric form such that its eigenvalue equations are much more easier to handle. 

i-i 



u 3 = Hu(2p-l,l;2p,l\2iK)]lu(2p,C;2p-l,C+l\2iK) 

I A B • • • -fiM 
fit A B ••• 
fit ^4 £ o ••• 

••• fit A B 
-B ••• fit -4 



(93) 



(94) 



2mlx2ml 



A 



B 



c / 



2mx2m 



o o ... 

is 



: : is 













V 



fit 



-is ••• 

••• 

-is 

; ; ; o o 



o o 

is 



2mx2m 

Now we proceed to solve the eigenvalue equation 



-is 




2mx2m 



(96) 



(97) 
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where u = {2sinh2K) ml / 2 uo z [uo\ /2 oo 2 ooi% 



m = 



1 z4> ^ 

z 2 i> 

Z 3 l/j Q 



( y u \ 

y 3 u 



u = 



V y u j 



Ml 
U 2 



By imposing the constraint, 

z l = -1 , orz = e'T (i 2 = 1,3,5- ••), 
one reduces eigenvalue equation (97) to 

(A + zB + z- 1 B ] ) n^j = \ij . 
Further, imposing the constraint, 

y m = -1, or y = (h = 1,3,5-..), 

(100) is reduced to 



V( A +y B + y~ 1 B r ) u = Xu, 



where 



V 



c —iz 1 s 



IZ s 



(98) 



(99) 

(100) 

(101) 
(102) 
(103) 



(104) 



c c — cos 



?s*c - ?(— Y^)se ™ -2(-^)se l ™ 



-is*c + i(l±^) se ^ + i(^±^)se 



-l+c*^-i- 
2 



C C — COS 



The eigenvalue equation (102) is further reduced to 

A 2 - 2 cosh 7 A + 1 = . (105) 

A = exp ( ± 7) is the solution of A . 7 is determined by 

c 3 

cosli7 = c(cos#i + cos#2) + sc cos#i cos^2 + s 2 sin#i sin0 2 , (106) 

5 



where two continuous variables, 9 1 , 9 2 , are obtained by taking the thermodynamic limit, 

m,Z->oo, 2£->0i, 2f*->0 2 . 
The partition function of V is 



(107) 
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The free energy per site under the thermodynamic limit is 

/ = -- lim N-HnZ (108) 

= -I In (2sinh2K) 1 / 2 - -L t 7 (9 1 , 9 2 )d9 1 d9 2 , (109) 

(J 6TT z p JO JO 

where (3 = j,. For convenience, J = 1, the internal energy per site is 

« = ^(Af) (no) 

1 /- 27r r 2 ^ 9t 

= -coth2K-— / / ^d9 1 d9 2 . (Ill) 
87H Jo io OA 

5 High Temperature Limit 

Let us expand u in terms of x at high temperature limit, x small, and compare the series 
of u obtained from the computer graphic method. @ 

1 + x 2 2x 

x= tanhAT, c= cosh2K = -, s= sinh2ir= -, (112) 

1 — x 2 1 — X 2 

(1+x 2 ) 3 l + x 2 , Q 
(l-x 2 )2x l-x 2 

■ COS &i COS 9 2 H 2 Sm ^l cos "2, (,11-j) 



(l-x 2 ) 2 1 2 (l-x 2 ) 2 
97 , , x 1 r (1 + x 2 ) 2 (l - 10x 2 + x 4 ) 2a; , „ 



2 J 



l + 6x 2 + x 4 . _ 4x(l + a; 2 ) ..... 
+ ^,^2)2 cosflicos0 2 + v _ y singising 2 }. (114) 

With the help of computer program Maple V Release 5.1, we get our result 

u = -3x- 8x 3 - 28x 5 - 132x 7 - 832x 9 + 0(x n ). (115) 

Comparing the result obtained by the computer graphic method, the series expansion of 
u for a simple cubic lattice can be transformed from the partition function, 

= 2cosh 3 K (1 + 3x 4 + 22x 6 + 187.5s 8 + 1980a; 10 + 0(x 12 )), (116) 



so 



u = -q^^Z* (H7) 



-3x - 12x 3 - 120x 5 - 1332x 7 - 17676x 9 + 0{x u ). (118) 
13 



(115 ) and (118) have very similar structures. The coefficient of x^ 1 does not exist in our 
result (115) though it exists naively in coth2K of (111). The first nonvanishing term — 3a; 
of (115) and the vanishing of all coefficients of even powers of x are the same as (118 ). All 
terms, up to 0(x 9 ), with minus signs in (115) are also the same as (118). The reason why 
(115) and (118) do not have the same first three or four coefficients may come from the 
fact that we cannot decipher precisely the largest eigenvalue of V . On the other hand, if 
we choose different directions for the 2D rotations, the eigenvalue equations may not be 
reduced to only one equation. This is a dilemma. 

The correct relation of cosh 7 implies the correct thermodynamic quantities of the 3D 
Ising model. To guess the correct relation as (106 ) may be a promising way for finding 
the right resolution of the 3D Ising model. For example, trying to repair (106), if we 
multiply the matrices A, B and B\ or a; 3 , with a correcting factor $(x), then we have 

c 3 

cosli7 = ccos#i + $(— ccos#2 + sc cos#i cos 62 + s 2 sin#i sm9 2 ). (H9) 

s 

With the help of Maple V, setting 

$( x ) = 1 _ ^ _ ^.4 _ 249^ _ 12325^ _ Q ^ w ^ ^ 

then the series expansion of u has the same result as (118). $(x) may be explained as, 
something like, a weighting function for different directions of 2D rotations. 
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